%=== This function solves for the model's value and policy functions
%=== It also solves for the steady state distributions of agents across
%=== possible states
function [A_staffing,VO,VN,U,U_hat,opt_ind_hp_N,opt_ind_hp_O,GU,GN,GO,opt_p_U,opt_p_N,opt_p_O,opt_x_U,opt_x_N,opt_x_O,opt_eta_U,opt_eta_N,opt_eta_O,opt_theta_U,opt_theta_N,opt_theta_O,dist_U,dist_N,dist_O] = DO_main_func(param,indicator_policy)

    %----- Specify Parameters -----%
    nu = param(1);
    beta = param(2);            % Discount factor (factors in exit probability)

    alpha = param(3);           % Production function parameter
    d = param(4);               % Human capital depreciation rate
    ell = param(5);             % Matching function parameter
    lambda = param(6);          % Probability of search OTJ
    chi = param(7);             % Investment cost parameter
    b=param(8);                 % Unemployment benefit

    muA =param(9);              % Staffing agency fee
    c =param(10);               % Cost of posting a vacancy
    delta =param(11);           % Separation rate
    h_sigma =param(12);         % Standard deviation of initial human capital distribution

    pr_phi(1) = param(13);      % Probability of entering the model with phi=1 (low fixed productivity type)
    pr_phi(2) = param(14);      % Probability of entering the model with phi=2 (high fixed productivity type)

    z(1) = param(15);           % Productivity parameter f(phi,h) = z(phi)h^(alpha)
    z(2) = param(16);

    muA_param=muA;
    gamma=2;
    vfi_toler = 0.0001;
    dist_toler = 0.0000001;
    
    % Specify grids %
    nphi = 2;                  % Number of possible phi values (fixed productivity types)
    hlb = 0.0;                 % Human capital grid lower bound 
    hub = 7.0;                 % Human capital grid upper bound 
    nh = 3000;                 % Number of possible human capital values
    h_grid = linspace(hlb,hub,nh); % Create human capital grid
    
    % ======================= Create discretized h distribution ==================== %
    h_mu = 1;                  % Mean of initial human capital distribution
    cdf_h = zeros(1,nh);       % CDF of initial h distribution
    pdf_h = zeros(1,nh);       % PDF of initial h distribution

    h_bin_size = (hub-hlb)/(nh-1);
    for ind_h = 1:nh
        if (ind_h ==1)
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
        else
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            bin_ub_last = h_grid(ind_h-1) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma) - normcdf(bin_ub_last,h_mu,h_sigma);
        end
    end

    %% ============== Value Function Iteration to Get Value and Policy Functions ================ %%
    %---- Prespecify Value Functions
    
    U = zeros(nphi,nh);                % Value of unemployment in production phase                     
    U_hat  = zeros(nphi,nh);           % Value of unemployment in search phase - when searching for optimal job type
    U_hat_update  = zeros(nphi,nh);
    
    VN = zeros(nphi,nh);               % Value of employment in non-outsourced (N) job in production phase
    VO = zeros(nphi,nh);               % Value of employment in outsourced (O) job in production phase
    VN_hat = zeros(nphi,nh);           % Value of employment in N job in search phase
    VO_hat = zeros(nphi,nh);           % Value of employment in O job in search phase
    VN_hat_update = zeros(nphi,nh);
    VO_hat_update = zeros(nphi,nh);
    
    A_staffing = zeros(nphi,nh);
    
    %---- Prespecify Value Functions
    opt_ind_hp_N = ones(nphi,nh);     % Optimal index for next period h value (h') if in non-outsourced job
    opt_ind_hp_O = ones(nphi,nh);     % Opitmal index for next period h value (h') if in outsourced job
    
    GU = zeros(nphi,nh);              % Indicates optimal job type to search for if unemployed - 1=non-outsourced (N) and 2=outsourced (O)
    GN = zeros(nphi,nh);              % Indicates optimal job type to search for if employed in N job - 1=non-outsourced and 2=outsourced
    GO = zeros(nphi,nh);              % Indicates optimal job type to search for if employed in O job - 1=non-outsourced and 2=outsourced
    
    opt_p_U = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when unemployed
    opt_p_N = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in N job
    opt_p_O = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in O job
    
    opt_x_U = zeros(nphi,nh);         % x contract value associated with optimal search choice when unemployed
    opt_x_N = zeros(nphi,nh);         % x contract value associated with optimal search choice when employed in N job
    opt_x_O = zeros(nphi,nh);         % x contract value associated with optimal search choice when employed in O job
    
    opt_eta_U = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when unemployed
    opt_eta_N = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in N job 
    opt_eta_O = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in O job
    
    opt_theta_U = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when unemployed
    opt_theta_N = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in N job 
    opt_theta_O = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in O job
    
    % Find index for next period human capital (h') if the agent only lets human capital depreciate
    ind_hp_min = ones(1,nh);          % Minimum index for h' (index for h' if human capital only depreciates)
    for ind_h = 1:nh
        h = h_grid(ind_h);
        error_min = 100;
        for ind_hp = 1:nh
            hp = h_grid(ind_hp);       % hp is the value of h'
            error = abs(hp - h*(1-d)); % Find index that minimizes (hp-h*(1-d))
            if (error<error_min)
                error_min = error;
                ind_hp_min(ind_h) = ind_hp;
            end
        end
    end
    
    %--------------- Begin Value Function Iteration ----------------%
    vfi_error = 1; % Set vfi_error > vfi_toler
    while (vfi_error > vfi_toler)
    
        for ind_phi = 1:nphi          % Loop over phi types
            for ind_h = 1:nh          % Loop over current h values
                h = h_grid(ind_h);    % h = current human capital (h) value 
                
                b_effect = b*z(ind_phi)*(h^alpha);  % Unemployment value
                c_effect = c;                                % Search cost
                
                %===== For section 5.4 =====% (all other policy changes only affect parameters stored in param)%

                if (indicator_policy==2)      % policy = 2: Outsourcing only in non-college jobs
                    if (ind_phi ==2)
                        muA = 1;
                    else
                        muA= muA_param;
                    end
                elseif (indicator_policy==3)   % policy = 3: Outsourcing only in college job
                    if (ind_phi ==1)
                        muA = 1;
                    else
                        muA= muA_param;
                    end 
                elseif (indicator_policy==4)  % policy = 4: Outsourcing only in h>=0.5
                    if (h >1)
                        muA = 1;
                    else
                        muA= muA_param;
                    end
                elseif (indicator_policy==5)  % policy = 5: Outsourcing only in h>=0.75
                    if (h >3.5)
                        muA = 1;
                    else
                        muA= muA_param;
                    end
                elseif (indicator_policy==6)  % policy = 6: Outsourcing only in h>=1
                    if (h >4)
                        muA = 1;
                    else
                        muA= muA_param;
                    end
                elseif (indicator_policy==7)  % policy = 7: Outsourcing only in h>=1.25
                    if (h >4.5)
                        muA = 1;
                    else
                        muA= muA_param;
                    end
                elseif (indicator_policy==8)  % policy = 8: Outsourcing only in h>=1.5
                    if (h >5)
                        muA = 1;
                    else
                        muA= muA_param;
                    end
                else
                end
                   
                % ---------- Solve hc Investment Decision ----------%
                max_VN = -9999;
                max_VO = -9999;
    
                for ind_hp = ind_hp_min(ind_h):nh      % Loop over possible next period hp (h') values
                    hp = h_grid(ind_hp);
    
                    temp_VN = z(ind_phi)*(h^alpha) - ((hp - h*(1-d))*chi)^gamma ...
                        + beta*(1-delta)*VN_hat(ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp);
    
                    temp_VO = (1-muA)*(z(ind_phi)*(h^alpha)) - ((hp - h*(1-d))*chi)^gamma ...
                        + beta*(1-delta)*VO_hat(ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp);
    
                    if (temp_VN>max_VN)                % Save ind_hp where N job employment value is highest
                        max_VN = temp_VN;
                        VN(ind_phi,ind_h) = temp_VN;
                        opt_ind_hp_N(ind_phi,ind_h) = ind_hp;
                    end
                    if (temp_VO>max_VO)                % Save ind_hp where O job employment value is highest
                        max_VO = temp_VO;
                        VO(ind_phi,ind_h) = temp_VO;
                        opt_ind_hp_O(ind_phi,ind_h) = ind_hp;
                    end
    
                end % end loop over hp options
                U(ind_phi,ind_h) = b_effect + beta*U_hat(ind_phi,ind_hp_min(ind_h));
                
                ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                A_staffing(ind_phi,ind_h) = muA*z(ind_phi)*(h^alpha) + (1-delta)*beta*(1-lambda*opt_p_O(ind_phi,ind_h))*A_staffing(ind_phi,ind_hp);

    
                %---------------------------------------------------%
                % ------- Solve Search Decision of Unemployed ----- %
                xO = VO(ind_phi,ind_h) - c_effect;
                etaO = (xO - U(ind_phi,ind_h))/(VO(ind_phi,ind_h) - U(ind_phi,ind_h));
                t_Uhat_O = xO;
    
                thetaN = (((VN(ind_phi,ind_h) - U(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                if (~isreal(thetaN))
                    thetaN = 0;
                    pN = 0;
                    xN = 0;
                    etaN = 0;
                else
                    pN = thetaN/((1+thetaN^ell)^(1/ell));
                    qN = pN/thetaN;
                    xN = VN(ind_phi,ind_h) - (c_effect/qN);
                    etaN = (xN - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                end
                t_Uhat_N = pN*xN + (1-pN)*U(ind_phi,ind_h);
    
                % Now Consider which job type's value is highest
                if (t_Uhat_N>t_Uhat_O)
                    GU(ind_phi,ind_h) = 1;
                    opt_p_U(ind_phi,ind_h) = pN;
                    U_hat_update(ind_phi,ind_h) = t_Uhat_N;
                    opt_x_U(ind_phi,ind_h) = xN;
                    opt_eta_U(ind_phi,ind_h) = etaN;
                    opt_theta_U(ind_phi,ind_h) = thetaN;
                elseif (t_Uhat_O>t_Uhat_N)
                    GU(ind_phi,ind_h) = 2;
                    opt_p_U(ind_phi,ind_h) = 1;
                    U_hat_update(ind_phi,ind_h) = t_Uhat_O;
                    opt_x_U(ind_phi,ind_h) = xO;
                    opt_eta_U(ind_phi,ind_h) = etaO;
                    opt_theta_U(ind_phi,ind_h) = 1.0;
                else
                    GU(ind_phi,ind_h) = 3;
                    opt_p_U(ind_phi,ind_h) = pN;
                    U_hat_update(ind_phi,ind_h) = t_Uhat_N;
                    opt_x_U(ind_phi,ind_h) = xN;
                    opt_eta_U(ind_phi,ind_h) = etaN;
                    opt_theta_U(ind_phi,ind_h) = thetaN;
                end
    
                %---------------------------------------------------%
                % ------- Solve Search Decision of Employed --------%
                xO = VO(ind_phi,ind_h) - c_effect;
                etaO = (xO - U(ind_phi,ind_h))/(VO(ind_phi,ind_h) - U(ind_phi,ind_h));
    
                t_VNhat_O = lambda*xO + (1-lambda)*VN(ind_phi,ind_h);
                t_VOhat_O = lambda*xO + (1-lambda)*VO(ind_phi,ind_h);
    
                thetaN_N = (((VN(ind_phi,ind_h) - VN(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                if (~isreal(thetaN_N))
                    thetaN_N = 0;
                    pN_N = 0;
                    xN_N = 0;
                    etaN_N = 0;
                else
                    pN_N = thetaN_N/((1+thetaN_N^ell)^(1/ell));
                    qN_N = pN_N/thetaN_N;
                    xN_N = VN(ind_phi,ind_h) - (c_effect/qN_N);
                    etaN_N = (xN_N - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                end
                t_VNhat_N = lambda*pN_N*xN_N + (1-lambda*pN_N)*VN(ind_phi,ind_h);
    
                thetaO_N = (((VN(ind_phi,ind_h) - VO(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                if (~isreal(thetaO_N))
                    thetaO_N=0;
                    pO_N = 0;
                    xO_N = 0;
                    etaO_N = 0;
                else
                    pO_N = thetaO_N/((1+thetaO_N^ell)^(1/ell));
                    qO_N = pO_N/thetaO_N;
                    xO_N = VN(ind_phi,ind_h) - (c_effect/qO_N);
                    etaO_N = (xO_N - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                end  
                t_VOhat_N = lambda*pO_N*xO_N + (1-lambda*pO_N)*VO(ind_phi,ind_h);
    
                % For those Currently in N - Now Consider which job type is highest
                if (t_VNhat_N>t_VNhat_O)
                    GN(ind_phi,ind_h) = 1;
                    opt_p_N(ind_phi,ind_h) = pN_N;
                    VN_hat_update(ind_phi,ind_h) = t_VNhat_N;
                    opt_x_N(ind_phi,ind_h) = xN_N;
                    opt_eta_N(ind_phi,ind_h) = etaN_N;
                    opt_theta_N(ind_phi,ind_h) = thetaN_N;
                elseif (t_VNhat_O>t_VNhat_N)
                    GN(ind_phi,ind_h) = 2;
                    opt_p_N(ind_phi,ind_h) = 1;
                    VN_hat_update(ind_phi,ind_h) = t_VNhat_O;
                    opt_x_N(ind_phi,ind_h) = xO;
                    opt_eta_N(ind_phi,ind_h) = etaO;
                    opt_theta_N(ind_phi,ind_h) = 1.0;
                else
                    GN(ind_phi,ind_h) = 3;
                    opt_p_N(ind_phi,ind_h) = pN_N;
                    VN_hat_update(ind_phi,ind_h) = t_VNhat_N;
                    opt_x_N(ind_phi,ind_h) = xN_N;
                    opt_eta_N(ind_phi,ind_h) = etaN_N;
                    opt_theta_N(ind_phi,ind_h) = thetaN_N;
                end
    
                % For those Currently in O - Now Consider which job type is highest
                if (t_VOhat_N>t_VOhat_O)
                    GO(ind_phi,ind_h) = 1;
                    opt_p_O(ind_phi,ind_h) = pO_N;
                    VO_hat_update(ind_phi,ind_h) = t_VOhat_N;
                    opt_x_O(ind_phi,ind_h) = xO_N;
                    opt_eta_O(ind_phi,ind_h) = etaO_N;
                    opt_theta_O(ind_phi,ind_h) = thetaO_N;
                elseif (t_VOhat_O>t_VOhat_N)
                    GO(ind_phi,ind_h) = 2;
                    opt_p_O(ind_phi,ind_h) = 1;
                    VO_hat_update(ind_phi,ind_h) = t_VOhat_O;
                    opt_x_O(ind_phi,ind_h) = xO;
                    opt_eta_O(ind_phi,ind_h) = etaO;
                    opt_theta_O(ind_phi,ind_h) = 1.0;
                else 
                    GO(ind_phi,ind_h) = 3;
                    opt_p_O(ind_phi,ind_h) = pO_N;
                    VO_hat_update(ind_phi,ind_h) = t_VOhat_N;
                    opt_x_O(ind_phi,ind_h) = xO_N;
                    opt_eta_O(ind_phi,ind_h) = etaO_N;
                    opt_theta_O(ind_phi,ind_h) = thetaO_N;
                end
    
            end % End loop over h values
        end % End loop over phi values
    
        error_U_hat = max(max(abs(U_hat_update(:,:) - U_hat(:,:))))/max(max(abs(U_hat_update(:,:))));
        error_VN_hat = max(max(abs(VN_hat_update(:,:) - VN_hat(:,:))))/max(max(abs(VN_hat_update(:,:))));
        error_VO_hat = max(max(abs(VO_hat_update(:,:) - VO_hat(:,:))))/max(max(abs(VO_hat_update(:,:))));
    
        vfi_error = max([error_U_hat,error_VN_hat,error_VO_hat]);
    
        
        U_hat(:,:) = U_hat_update(:,:);
        VN_hat(:,:) = VN_hat_update(:,:);
        VO_hat(:,:) = VO_hat_update(:,:);
        
    end % End while (vfi_error > vfi_toler)  
    
    %==================================================================================================================%
    %% ================================ Solve for Steady State Distribution ========================================= %%
    %==================================================================================================================%
    
    dist_U = zeros(nphi,nh);           % Steady state mass of unemployed agents
    dist_N = zeros(nphi,nh);           % Steady state mass of agents employed in non-outsourced job
    dist_O = zeros(nphi,nh);           % Steady state mass of agents employed in outsourced job
    
    dist_U_update = zeros(nphi,nh);
    dist_N_update = zeros(nphi,nh);
    dist_O_update = zeros(nphi,nh);
    
    dist_U(1,1) = pr_phi(1);           % Start with a guess regarding how the mass of agents is distributed
    dist_U(2,1) = pr_phi(2);
    
    dist_error = 1;
    while (dist_error>dist_toler)      % Iterate until distribution guess equals resulting distribution
    
        for ind_phi = 1:nphi           % Loop over phi types
            for ind_h = 1:nh           % Loop over current h values 
                %=================================%
                %========= New Entrants ==========%
                %=================================%
                % New Entrant - Finds Job
                if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)     % Searches for N
                    dist_N_update(ind_phi,ind_h) = dist_N_update(ind_phi,ind_h)...
                        + opt_p_U(ind_phi,ind_h)*pr_phi(ind_phi)*nu*pdf_h(ind_h);
                else                                                   % Searches for O
                    dist_O_update(ind_phi,ind_h) = dist_O_update(ind_phi,ind_h)...
                        + opt_p_U(ind_phi,ind_h)*pr_phi(ind_phi)*nu*pdf_h(ind_h);
                end
                % New Entrant - Does Not Find Job
                dist_U_update(ind_phi,ind_h) = dist_U_update(ind_phi,ind_h) + (1-opt_p_U(ind_phi,ind_h))*pr_phi(ind_phi)*nu*pdf_h(ind_h);
    
                %=================================%
                %===== Previously Unemployed =====%
                %=================================%
                ind_hp = ind_hp_min(ind_h);             % New h value for the unemployed is what they get from allowing human capital to depreciate
    
                % Unemployed - Finds a job
                if (GU(ind_phi,ind_hp) == 1 || GU(ind_phi,ind_hp) == 3)     % Searches for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*dist_U(ind_phi,ind_h)*(1-nu);
                else                                                         % Searches for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*dist_U(ind_phi,ind_h)*(1-nu);
                end 
                % Unemployed - Does not find a job
                dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp) + (1-opt_p_U(ind_phi,ind_hp))*dist_U(ind_phi,ind_h)*(1-nu);
    
                %=================================%
                %======= Previously in N =========%
                %=================================%
                ind_hp = opt_ind_hp_N(ind_phi,ind_h);   % New h value for employed in N job is the result of their optimal investment choice
    
                % ----- N - No Separation Shock - Finds New Job
                if (GN(ind_phi,ind_hp)==1 || GN(ind_phi,ind_hp)==3)      % N - Search OTJ for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + lambda*opt_p_N(ind_phi,ind_hp)*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
                else                                                      % N - Search OTJ for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + lambda*opt_p_N(ind_phi,ind_hp)*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
                end
    
                % ----- N - No Separation Shock - Does Not Find New Job
                dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                    + (1-lambda*opt_p_N(ind_phi,ind_hp))*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
    
                % ----- N - Separation Shock - Finds Job
                if (GU(ind_phi,ind_hp)==1 || GU(ind_phi,ind_hp)==3)       % N - Separation - Search for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*delta*dist_N(ind_phi,ind_h)*(1-nu);
                else                                                       % N - Separation - Search for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*delta*dist_N(ind_phi,ind_h)*(1-nu);
                end
                
                % ----- N - Separation Shock - Does Not Find New Job
                dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp)...
                    + (1-opt_p_U(ind_phi,ind_hp))*delta*dist_N(ind_phi,ind_h)*(1-nu);
    
                %=================================%
                %======= Previously in O =========%
                %=================================%
                ind_hp = opt_ind_hp_O(ind_phi,ind_h);   % New h value for employed in O job is the result of their optimal investment choice
    
                % ----- O - No Separation Shock - Finds New Job
                if (GO(ind_phi,ind_hp)==1 || GO(ind_phi,ind_hp)==3)      % O - Search OTJ for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + lambda*opt_p_O(ind_phi,ind_hp)*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
                else                                                      % O - Search OTJ for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + lambda*opt_p_O(ind_phi,ind_hp)*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
                end
                % ----- O - No Separation Shock - Does Not Find New Job
                dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                    + (1-lambda*opt_p_O(ind_phi,ind_hp))*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
    
                % ----- O - Separation Shock - Finds Job
                if (GU(ind_phi,ind_hp)==1 || GU(ind_phi,ind_hp)==3)       % O - Separation - Search for N
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*delta*dist_O(ind_phi,ind_h)*(1-nu);
                else                                                       % O - Separation - Search for O
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + opt_p_U(ind_phi,ind_hp)*delta*dist_O(ind_phi,ind_h)*(1-nu);
                end
    
                % ----- O - Separation Shock - Does Not Find New Job
                dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp)...
                    + (1-opt_p_U(ind_phi,ind_hp))*delta*dist_O(ind_phi,ind_h)*(1-nu);
    
            end % End loop over last period h value
        end % End loop over phi values
    
        error_dist_U = max(max(abs(dist_U_update(:,:) - dist_U(:,:))));
        error_dist_N = max(max(abs(dist_N_update(:,:) - dist_N(:,:))));
        error_dist_O = max(max(abs(dist_O_update(:,:) - dist_O(:,:))));
    
        dist_error = max([error_dist_U,error_dist_N,error_dist_O]);
    
        dist_U(:,:) = dist_U_update(:,:);
        dist_N(:,:) = dist_N_update(:,:);
        dist_O(:,:) = dist_O_update(:,:);
    
        dist_U_update(:,:) = 0;
        dist_N_update(:,:) = 0;
        dist_O_update(:,:) = 0;
    
    
    end % end while (dist_error>dist_toler)
    
    %dist_sum = sum(sum(dist_U(:,:))) + sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:)));


end % End function